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Numerical simulations are performed of the approach to the singularity in Gowdy spacetimes on 
S^ X S^ X R. The behavior is similar to that of Gowdy spacetimes on T"^ x R. In particular, the 
singularity is asymptotically velocity term dominated, except at isolated points where spiky features 
develop. 
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Q\ ' I. INTRODUCTION 

a^■ 

There have been several numerical investigations of the approach to the singularity in inhomogeneous cosmologies. 
-|[ In general, it is found that (except at isolated points) the approach to the singularity is either asymptotically 
C ' velocity term dominated |^ (AVTD) or is oscillatory. In the oscillatory case, there are epochs of velocity term 
dominance punctuated by short "bounces." The most extensively studied inhomogeneous cosmology is the Gowdy 
spacetime ||l^,|ll| on T^ x R. Here the approach to the singularity is AVTD except at isolated points. The Gowdy 
spacetimes on T'^ x _R are especially well suited to a numerical treatment for the following reasons: (i) Due to the 
presence of two Killing fields, the metric components depend on only two spacetime coordinates, (ii) The constraint 
^ . equations are easy to implement, (iii) The boundary conditions are particularly simple, just periodic boundary 
^\^ • conditions in the one nontrivial spatial direction. 

The original work of Gowdy |ll[| treated spatially compact spacetimes with a two parameter spacelike isometry 
^D ' group. Gowdy showed that, for these spacetimes, the topology of space must be T^ or S^ or S^ x S^. Given the 
>0 numerical results for the T^ case, it is natural to ask what happens in the other two cases. In a recent paper |lj] 
Obregon and Ryan note that the Kerr metric between the outer and inner horizons is a Gowdy spacetime with spatial 
topology S^ X S^. They analyze the behavior of this spacetime and speculate that there may be significant diff'erences 
between the behavior of Gowdy spacetimes on T'^ x i? and on S^ x S^ x R. 
y^ ,' A numerical simulation of the S^ x S^ case presents some difficulties that are absent in the T"^ case. The constraint 

equations become more complicated, and there are difficulties associated with boundary conditions. In the T^ case. 
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^ ■ the Killing fields are nowhere vanishing. However, in the 5^ x S^ case, one of the Killing fields vanishes at the north 
rrj and south poles of the S^. Smoothness of the metric at these axis points then requires that the metric components 
I> . behave in a particular way at these points. A computer code to evolve the S'^ x S^ case therefore must enforce these 
t.^ ' smoothness conditions as boundary conditions, and must do so in such a way that the evolution is both stable and 
'^j , accurate. These issues are similar to those encountered in the numerical evolution of axisymmetric spacetimes, and 
d ' the techniques presented here for Gowdy spacetimes should be useful for axisymmetric spacetimes as well. 

This paper presents the results of numerical simulations of Gowdy spacetimes on S^ x S^ x R. Section 2 presents 
the metric and vacuum Einstein equations in a form suitable for numerical evolution. The numerical technique is 
presented in section 3, with the results given in section 4. 

II. METRIC AND FIELD EQUATIONS 

The Gowdy metric on 5*^ x S"^ x i? has the form ||ri| 

ds"^ = e*^ {-dt^ + d9^) + sini sin 6* (e^[d(j) + Q dSf + e'^dS^') . (1) 

Here the metric functions M, L and Q depend only on t and 6. Thus our two Killing fields are {d/d(f) and {d/dS). 
The coordinates cj) and S are identified with period 27r, with S the coordinate on the S^ and {0, <f>) the coordinates on 
the 5*2. 
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Here the "axis" points are at 6* = and 9 = tt. The spacetime singularities ("big bang" and "big crunch") are 
at i = and t = tt. This form of the metric presents difficulties for a numerical treatment. Smoothness at the axis 
requires divergent behavior in the functions L and M. Furthermore, the spactimes singularities occur at finite values 
of the time coordinate. This is likely to lead to bad behavior of the numerical simulation near t = or i = tt. These 
difficulties are overcome with a new choice of metric functions and time coordinate. Define the new metric functions 
P and 7 by 



P = L - lnsin6' 



(2) 



Define the new time coordinate r by 



The metric then takes the form 
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Smoothness of the metric at the axis is equivalent to the requirement that P, Q and 7 be smooth functions of cos 9 
with 7 vanishing at 6* = and 9 = tt. Note that for / any smooth function of cos 6', it follows that df/d9 = at 6* = 
and 9 = TT. 

As in the T^ case, the vacuum Einstein field equations become evolution equations for P and Q and "constraint" 
equations that determine 7. The evolution equations are 
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Here a subscript denotes partial derivative with respect to the corresponding coordinate. Note that, as in the T^ case, 
the evolution equations have no dependence on 7. 
The constraint equations are 



cot 9^T — tanh t ^e = A 



(8) 
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cosh r 
where the quantities A and B are given by 

2^EEtanhTPe + Pr Pe + e^^ 8111^9 QrQe 



AB ~ 2 t&iihT Pr + [Prf + e^^ sin26' (g^)^ + tanhV - 4 



Solving equations (ph and (0) for 79 and 7^- we find 
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Given a solution of the evolution equations (0) and (M) for P and Q, equations ( 12) a nd (O) and the smoothness 
condition that 7 = at = completely determine 7. Actually, equations (|lj) and (|l3|) seem to be in danger of over 



determining 7, but the integrability condition for these equations is automatically satisfied as a consequence of the 
evolution equations for P and Q. There are, however, two remaining difficulties with the equations for 7. The first 
has to do with the fact that the denominator in equations (^^ and (13) vanishes when | cot6\ — sinhr. Smoothness 
of the metric then requires that the numerators of these equations vanish whenever the denominator does. This places 
conditions on P and Q. If these conditions are satisfied for the initial data, the evolution equations will preserve them. 
The second difficulty has to do with the fact that 7 must vanish at 6 ~ n as well as 6 = 0. Integrating equation (flT 
from to vr it then follows that wc must have 



cosher (A tanhr + B cot 9) d9 
cot'^9 — sinh^T 







(14) 



If this condition is satisfied by the initial data, then the evolution equations will preserve it. 

In summary, the initial data for P and Q are not completely freely specifiable. They must satisfy conditions at 
the points where | cot^j = sinhr as well as an integral condition. Given initial data satisfying these conditions, the 
evolution equations (o) and (R) then determine P and Q and the constraint equations (O) and (0) then determine 
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III. NUMERICAL METHODS 



We now turn to the numerical methods used to implement the evolution equations. We begin by casting the 
equations in first order form by introducing the quantities V = Pr and W = Qr- These quantities then satisfy the 
equations 



Vr = e^^ sin^ 9W^ 
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Thus the evolution equations have the form X^ = F{X,t). We implement these equations using an iterative 
Crank-Nicholson scheme. Given X at time r, we define Xq{t + At) = X{t) and then iterate the equation 



X„+i(r + Ar)=X(r) 



At 
~2~ 



F{X{t),t) + F(X„(T + Ar),T + Ar) 



(17) 



In principle, one should iterate until some sort of convergence is achieved. In practice, we simply iterate 10 times. 
We use At — A9/2 where /S.9 is the spatial grid spacing. 

The spatial grid is as follows: Let ng be the number of spatial grid points. Then we choose /S.9 — ^/{ng — 2) 
and 9i = {i — 1.5)A6'. Thus in addition to the "physical zones" for i = 2, 3, . . . , ne — 1, we have two "ghost zones" 
at 9i = — A0/2 and 9ne = tt + (A6'/2). The ghost zones are not part of the spacetime: variables there are set by 
boundary conditions. For any quantity S*, define Si = S{9i). Spatial derivatives are implemented using the usual 
second order scheme: 
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Smoothness of the metric requires that Pg = at 6* = 0. Since 9 = is halfway between i = 1 and i = 2, we 
implement this condition as Pi = P2. Similarly, we use Qi — Q2 since Qg = at = 0. Correspondingly, the 
requirement that Pg and Qg vanish at ^ = tt is implemented as Pn^ = Png-i and Q„g = Qn^-i- These boundary 
conditions are imposed at each iteration of the Crank-Nicholson scheme. 



IV. RESULTS 

To test the eomputer code, it is helpful to have some closed form exact solution of the evolution equations to 
compare to the numerical evolution of the corresponding initial data. In particular, for a second order accurate 
evolution scheme, the difference between the numerical solution and the exact solution should converge to zero as the 
grid spacing squared. 

A polarized Gowdy spacetime is one for which the Killing vectors are hypersurface orthogonal. For our form of the 
metric, that is equivalent to the condition Q = 0. For polarized Gowdy spacetimes, the evolution equation (R) for Q 
is trivially satisfied, and the evolution equation (^) for P reduces to the following: 

P^^ = — -^ [Peg + cotePe - 1] . (20) 

cosh T 

This is a linear equation which can be solved by separation of variables, though one must choose only those solutions 
that satisfy the additional conditions for smoothness of the metric. 

Unfortunately, the polarized solutions do not provide, by themselves, a very stringent code test: the evolution 
equation for Q and the nonlinear terms in the evolution equation for P are not tested at all. Fortunately, there is a 
technique, the Ehlers solution generating technique, which allows us to begin with a polarized solution and produce 
an unpolarized solution. Let P be any solution of the polarized equation ( |20| ) and let c be any constant. Define P 
and Q by 



P = P - In 
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— 2c 
Qr ^ — (2 cos 61 + sin6iPe) , (22) 
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Qe = 2csine' (tanhr - Pr) . (23) 

Then (P, Q) is a solution of the unpolarized Gowdy equations (g) and (0)- 
We use the following polarized solution: 

P = - Incoshr + 2t , (24) 

where 5 is a constant. The solution generating technique then yields an unpolarized solution with P given by equation 
( pl| ) and Q given by 

Q = 4c(l - tanhr) cose* . (25) 

Figure 1 shows P for the exact solution and the numerical evolution. (Here there are 502 spatial grid points, c = 1, 
and the initial data at r = are evolved to r = 10. The results are shown at 51 equally spaced points from 9 = to 
9 = it). Figure 2 shows the difference between exact and numerical solutions. Here the parameters are as in figure (1), 
except that two simulations are run: one with 502 spatial gridpoints and one with 1002 grid points. For comparison, 
the results on the finer grid are multiplied by a factor of 4. The results show second order convergence. (Note: due to 
the presence of the ghost zones, quantities must be interpolated on the grids to make a comparison between quantities 
at the same values of 9) . 

We would now like to find the generic behavior of Gowdy spacetimes on S*^ x S*^ x R. In the T^ case Q a family of 
initial data was chosen and evolved. It was argued that the behavior of these spacetimes reflects the generic behavior. 
Here, we choose a similar family. The initial data at r = are P = 0, Pr = vqCOS^, Q = 2cos^, Qr = 0. Here vo 
is a constant. These data satisfy the constraint conditions. Figures 3-8 show the evolution of these data for various 
values of the parameter vq. Here, wq = 2 in figures 3 and 4, wq = 4 in figure 5 and 6, and uq = 8 in figures 7 and 8. In 
all cases, the range of 9 is (0, tt), the range of r is (0, 10) and the simulation is run with 1002 spatial grid points. Note 
the presence of spiky features. The large r behavior of the solutions is the following: There are functions Qoo(^) and 
foo(^) with Vooi9) < 1 such that away from the spiky features we have Q — > Qoo{9) and Pr -^ foo(^) for large r. 

The reason for this behavior is not hard to find and is essentially the same as in the T^ case. For large t and 
provided that P is growing no faster than r, The terms in equations (0) and (R) proportional to 1/cosh r become 
negligible. The truncated equations obtained by neglecting these terms are 



P^^ ^e^^sin^ 0{Qr) , (26) 

Q,,^-2PrQr ■ (27) 

Equations (|2^) and ^^ are called the AVTD equations. They can be solved in closed form and have the property 
that Q -^ Qoo{&) and Pr -^ Voc{9) as r ^ oo. Solutions of the full equations (O) and (R) are called AVTD provided 
that they approach solutions of the AVTD equations for large r. Thus we have an explanation of the AVTD behavior 
provided that we can show that P grows no faster than r. As in the T^ case, if P grows faster than r then the term 
in equation (^) proportional to e^^'^/cosh r will cause a "bounce" that leaves P growing less fast than r. Thus, an 
analysis of the large r behavior of the evolution equations (H) and (Q), essentially the same as in the T^ case, leads to 
an explanation of the AVTD behavior. 

The AVTD behavior can also be explained by an analysis of the local properties of Gowdy spacetimes on S"^ x 5*^ x R. 
Define Sa = Va(sini sin0). Then in regions of the spacetime where Sa is timelike, the region is locally isometric to 
a Gowdy spacetime on T^ x R. In regions where Sa is spacelike, the region is locally isometric to a cylindrical wave. 
Thus, the behavior that we should expect in the S^ x S^ case is a combination of the the behavior of the T^ case 
and the behavior of cylindrical waves. Furthermore, for any point on the S^ except the poles, as the singularity is 
approached, Sa becomes timelike at that point. Thus the asymptotic behavior as the singularity is approached in the 
S"^ X S^ case should be the same as in the T^ case. 

We now turn to an analysis of the spiky features seen in the metric functions P and Q. The argument of the 
previous paragraph indicates that these features are essentially the same as those seen in the T^ case. In fact, these 
features can be explained using the evolution equations (^ and (^ as was done in the T^ case. Q For large r, it 
follows from equation (^ that Qr ~ IIq(^) e~^^ for some function IIq{9). Then, using this result in equation (^ we 
have an approximate evolution equation for P: 



Prr « sm^e 



{^qY r^ {Qef 
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These terms eventually drive Pr to the range between and 1. However, at a point 6i where Qg vanishes, Pt can be 
greater than 1. This leads to a spiky feature in P, since P^. > 1 at 6i but Pt- < 1 at points near 9i. This sort of spiky 
feature is illustrated in figure 9. Correspondingly, at a point 62 where Ilg vanishes, Pr can be less than zero. This 
leads to sharp features in P since Pr > at points near 02- Also since the region where P < leads to rapid growth 
in Q, there is a sharp feature in Q. This sort of feature is illustrated in figure 10. 

In summary, a numerical treatment of Gowdy spacetimes on 5^ x 5^ x P reveals that they are very similar to 
Gowdy spacetimes on T^ x R. In particular, they show the same behavior of AVTD behavior almost everywhere, and 
they have the same sort of spiky features at isolated points. 
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FIG. 1. The analytic values (line) and numerical values (dots) of P are plotted vs. 6 at r = 10. 
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FIG. 2. The diflferences between analytic and numerical values for P are compared for two different resolutions. The line is 
AP for 502 grid points and the squares are 4 times AP for 1002 grid points. The fact that these curves agree shows second 
order convergence 




FIG. 3. The evolution of P for initial data with vq = 2. 
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FIG. 4. The evolution of Q for initial data with vo = 2. 




FIG. 5. The evolution of P for initial data with vq = 4. 
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FIG. 6. The evolution of Q for initial data with vo — 4. 
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FIG. 7. The evolution of P for initial data with vq 
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FIG. 8. The evolution of Q for initial data with vo 
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FIG. 9. a spiky feature in P (solid line) occurs where Q (dashed line) has an extremum. Here, vo = 
simulation is run with 2002 spatial grid points. 
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FIG. 10. a sharp feature in Q (dashed Une) occurs at a downward spike in P (soUd hne). Here, wo = 4, r = 10 and the 
simulation is run with 2002 spatial grid points. 
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